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1. Introduction 

Symmetries play a crucial role in our understanding of strong interactions physics. Of par- 
ticular interest is the SU (Nf) x SU (Nf) chiral symmetry of the massless QCD Lagrangian. Its 
spontaneous and explicit symmetry breaking pattern is the starting point for a large part of low- 
energy phenomenology. Furthermore, an index theorem holds for the continuum Dirac operator 
which is believed to have significant consequences for low-energy physics too. This theorem states 
that the index v of the Dirac operator is equal to the topological charge Q top of the gauge field 

1 

v = n L -n R = -try$D(0) = Q ta p (1.1) 

with riL and hr respectively the number of left and right handed zero-modes of the massless Dirac 
operator D(0). Since these properties are important in the continuum, it is only natural to respect 
chiral symmetry when formulating QCD on the lattice. 

The correct form of chiral symmetry at finite lattice spacing [|l|] is the Ginsparg-Wilson (GW) 
relation [0] 

D(0)y 5 + 75 D(0) = -^-D(0)y 5 D(0) (1.2) 

with a the lattice spacing and Ro a parameter. For operators which fulfill the GW equation the 
index theorem ^ works at finite lattice spacing just like in the continuum. Moreover, the theory 
is automatically 0(a) improved, additive mass renormalization is absent and renormalization pat- 
terns can be greatly simplified [Qj. The fermion determinant is strictly positive if the quark mass 
is greater than zero. Prominent solutions of the Ginsparg-Wilson equation include domain wall 
fermions at large extent of the fifth dimension 01, the overlap operator [|7|] and the fixed point 
Dirac operator Unfortunately, the application of these operators is significantly more expen- 

sive than standard operators like the Wilson or staggered Dirac operator. This prevented so far a 
more wide-spread use. 

In recent years simulations of QCD on large lattices and at small quark masses have become 
possible. This is not only due to increased computer resources but rather because of algorithmic 



advances, see the reviews given at this conference[|10|, |1 These improvements apply to the sim- 
ulation of dynamical chiral fermions too. In one respect, however, chiral fermions are different. 
Because of the index theorem Eq. |1.1| , any chiral Dirac operator changes discontinuously where 
its index changes. This has to be accounted for in the equations of motion which are at the heart 
of most current algorithms. The different approaches how to deal with the change of topological 
sector algorithmically are the main subject of this write-up. 

This proceedings contribution is about dynamical simulations using Neuberger's overlap Dirac 
operator. However, also domain wall fermions at finite time extent, the parameterized fixed point 



operator [12] and chirally improved fermions [13], which all implement chiral symmetry to a high 



degree, have been used in dynamical simulations. Their approximate symmetry leads to an ap- 
proximation of the discontinuity encountered for actions with exact chiral symmetry. As we are 
going to discuss, this can cause large forces in the molecular dynamics evolution while changing 
topological sector. The issue of changing topology has therefore received considerable attention in 
all those simulations too. 
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This paper is organized as follows. First we review the definition of the overlap operator and 
the fermionic action with a focus on the discontinuity where the index changes. Three methods to 
deal with the discontinuity in the context of Hybrid Monte Carlo are discussed: modified molec- 



ular dynamics evolution in Sec. \^2\, approximation of the discontinuity in Sec. g3| and topology 



conserving actions in Sec. . The write-up closes with a summary. 
2. The action 

Let us start by reviewing the action whose simulation is the subject of this paper. The overlap 
Dirac operator [{7|] is given by 

D ov (m) = (R Q --)[l + y 5 e{h{-Ro))] +m (2.1) 

with m the bare quark mass. It is constructed from a Hermitian lattice Dirac operator h = y^d which 
is typically the standard Wilson operator. It is taken at the negative mass — Rq at the cut-off scale. 



Ro can be tuned to achieve optimal locality of the resulting overlap operator []14[]. In the following 
h(—Ro) is referred to as the kernel operator. e(h) is the matrix sign function. For Hermitian, 
non-singular matrices it can be defined by 

e(h) = ^= = ^sign(l i ) ¥i¥ J (2.2) 

with the sum running over all eigenvalues A, of the kernel operator h and Yi the associated eigen- 
modes. The numerical construction of e(h) usually uses a combination of the two definitions in 



Eq. 12. The lowest modes of h(—Ro) are computed explicitly and the spectral representation can 
be used. For the rest of the spectrum a polynomial or rational approximation to 1/ s/J? implements 
the sign function to very high accuracy. 

For simulations of two flavors one is also interested in H = D^D, with H = 75D. (Here 
and in the following we denote the overlap Dirac operator by capital letters D and H whereas the 
Hermitian kernel operator is denoted by small h.) Because of the Ginsparg-Wilson equation and 
75-Hermiticity, H 2 commutes with y 5 and is therefore block-diagonal in the chiralities with 

2 

H 2 a = P a H 2 P a = 2(R 2 - —)P a [1 + ae(h(-R Q ))]Pa + m 2 P a (2.3) 



where P a = ^(1 + CJ75) projects to chirality a = ±1. 

Obviously, the overlap operator is discontinuous where one of the eigenvalues Xq of the kernel 
operator changes sign. This is precisely where the topological charge as defined by the index of the 



Dirac operator changes. Evaluating the lattice version of Eq. 1.1 for the overlap operator [Bp gives 



V = i Tr7sD(m = 0) = ^ Tre (K-*o)) = ^EsignA; (2.4) 

with the sum running over all eigenvalues A; of the kernel operator. The direction of the sign change 
determines whether the topological charge decreases or increases by one unit. 

On the surfaces which separate the topological sectors, the spectrum of the overlap operator 
changes: One zero-mode is created or vanishes and also the rest of the spectrum moves. Interest- 
ingly, the ratio of the fermion determinants on the two sides of the step can be computed without 
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the explicit computation of the determinant itself: Because the eigenmodes of the kernel are con- 
tinuous functions of the gauge fields, the change of the sign of one of the eigenvalues amounts to 
the following change in the overlap operator 

D(m) =D{m)- (2R - m) sign(Ao) 75 VoVq ( 2 -5) 

with signAo the sign of the eigenvalue before the crossing and y/b the associated eigenfunction. 
This translates to a discontinuity of the fermion determinant JT^ ] 



detD(ra) 



1 - (2R - m) sign(Ao) yj,) 



detD(m) (2.6) 



which can be computed by just one inversion of the overlap operator on the crossing mode. How- 
ever, most of the time we do not deal with the determinant directly but introduce it into the func- 
tional integral by pseudo-fermion fields 

det// 2 oc j d0t e -^n-H = Jd$ d0 f e - s ^ u ^ (2.7) 

with Sf[$,U] the effective fermion action. The discontinuity in Sf[(j),U] can be computed using 
the Sherman-Morrison formula 



AS f = AyV-U= V asign(Ao) 1 — ^ — — i^r 1 2 

f k H l ahl 'l-asign(Ao)(4/?2-m2)v/ // ff 2 (m)^o IV H ff (m) 2 TO 

(2.8) 

which again can be evaluated by inverting the (squared) overlap operator on the crossing mode of 
the kernel operator. 

In summary, the overlap operator is constructed via the matrix sign function of a doubler free 
lattice Dirac operator. This makes its application roughly 20-200 times more expensive than a stan- 
dard Wilson operator. 1 The fermion determinant and the effective action for the overlap operator 
are non-zero and finite everywhere. However, it has discontinuities where the topological charge 
(the index) changes. They have their origin in one eigenvalue of the kernel operator changing sign. 
The height of these steps can be computed by one inversion of D^D on the crossing mode. 



3. Hybrid Monte Carlo 

The most popular exact algorithm used in dynamical fermion simulations is Hybrid Monte 
Carlo [17]. It is a combination of Molecular Dynamics (MD) Jl8| ] with a Metropolis accept/reject 



step that makes it exact. Like in all MD based simulations conjugate momenta n are introduced 
which drive the evolution of the generalized coordinates, the gauge fields in our case. The Hamilto- 
nian H[U,n,(j)] from which the equations of motion are derived is composed of the kinetic term for 
the momenta and the potential given by the action of the theory which we actually want to simulate 



H[U,n,<j>] = —+S f [U,<l>]+S g [U]. (3.1) 



'A recent comparison!^] finds the cost of computing the propagator with the overlap operator to be 30-120 times 
larger than with twisted mass fermions. 
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Sf[U,(j)] is the effective fermion action at fixed pseudo-fermion fields and S g [U] the gauge action. 

The elementary update step in the HMC algorithm is called a trajectory. At the beginning, 
a heat bath for the pseudo-fermion field and the momenta is performed. Then the equations of 
motion, which derive from dH /dl = are solved numerically for some (fictitious) time T. The 
simplest integrator is the leap frog which repeats 

[V — >T u (8r/2)U , n — >T p {8z)n , U — >T v {8z/2)U] (3.2) 

t/5t times; Ty and T p are the gauge field and momentum updates. At the end of the trajectory, an 
accept/reject step is performed which corrects for errors introduced by the numerical integration. 
The new gauge configuration is accepted with probability P acc = min{exp(— AH[U, 71, (j)}), 1}. 

The leap-frog integrator provides us with an approximate solution to the equations of motion 
with an error that vanishes with (St) 2 as 8x — ► 0. Decreasing the step-size can therefore bring 
the acceptance rate arbitrarily close to 1. This is not the case for the overlap action because the 
leap-frog (almost) never hits the discontinuity and its effect is therefore not taken into account in 
the numerical integration. Without a special treatment of the step in the action, the acceptance rate 
is therefore smaller than 1 ; typically far too small for any practical purpose. How to deal with this 
discontinuity is the subject of the rest of this write-up. In Sec. we discuss the modification of 
the MD evolution proposed by Fodor, Katz and Szabo (FKS) Jl9| , 20]. The subject of Sec. 3_3 are 
algorithms which approximate the step in the action and thereby make the evolution accessible to 
standard methods. This approximation has the to be corrected for to get the exact overlap action. 
Topology conserving actions, which solve the problem by avoiding the discontinuity completely, 



are described in Sec. 3.4. First, however, we discuss various possibilities to introduce the fermion 



determinant into the simulation. 



3.1 Multiple pseudo fermions 



In recent years, the quality of the estimator used to introduce the fermion determinant into 
the functional integral has turned out to be of crucial importance to the performance of standard 



algorithms. The pseudo-fermion field in Eq. 2.7 is such a stochastic estimator for the change in 



the fermion determinant over an update: At the starting configuration a pseudo-fermion heat-bath 
is performed, <p = D(m)R with R a Gaussian random field. The change in the determinant can be 
computed by the average over the random field 

det ^ 2 = / e -*Sf[M) (3 3) 

Using just one pseudo-fermion field for this estimate is known to be very noisy [^l|] which in the 
context of Hybrid Monte Carlo can cause large forces. 

A successful way of improvement is to replace the original fermion matrix by a product of 
matrices 

detg = detD j D = detMi detM 2 • --detM N . (3.4) 

Each of these determinants is then introduced by a pseudo-fermion field. If the matrices M, are 
suitably chosen, it can be easier to estimate the (change in) their determinant stochastically and the 
fluctuations are greatly reduced. 
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The initial idea by Hasenbusch was to use N = 2, M\ = Q[m)/Q(M) and M 2 = Q(M) with 



Q(M) the fermion matrix at a larger quark mass [ |22| , |23[ ] . Another choice is to use the A^-th root 



of Q for all M, ||24j, |25j]. So far neither of these two methods has proved superior to the other. The 
former, however, is easier to implement. As we will see in the following, the improved estimator 
not only reduces the fermionic forces but it can also drastically reduce the auto-correlation time of 



the topological charge [26]. 



3.2 Modified evolution 

Most of the experience with dynamical overlap fermions is based on the algorithm proposed 
by Fodor, Katz and Szabo (FKS) JT9[], which is used — sometimes in variants — in Refs. [26], 27, 28 
15 



, 30, 31 



, |33J]. So far the simulations have been limited to small volume, lattices up to 10 
sites albeit at a coarse lattice spacing. 

The FKS algorithm is the standard Hybrid Monte Carlo unless a discontinuity is encountered. 
This can only happen — an eigenvalue of the kernel operator h(—Ro) changes sign — during the 
gauge field update Tjj{8x/2). The corresponding update step is the split in two: one where the 
trajectory is advanced right to the point where the index changes by Tu(8%{), i.e. where the eigen- 
value of the kernel operator gets zero. There a discrete update of the momenta is performed. It 
depends on whether there is enough momentum 7i± perpendicular to the surface which separates 
the two topological sectors to get across the step in the action AS. In this case, the momentum is 
reduced such that energy is (approximately) conserved and the topological sector changes. If there 
is not enough momentum, n± is reversed, the trajectory is reflected off the surface and stays in 
the same topological sector. To be precise, the momentum component n± is changed to n' ± in the 
following way 



■ic± 



IAS 

[w] P 



if 2AS> |7Ti| 2 
if 2AS< \k ± \ 2 



(3.5) 



After that, the gauge fields are advanced to the end of the update step by 7j/(5t/2 — 8x\). This 
update is obviously reversible and was proved to be area conserving too Jl9|]. Combined with the 
standard Metropolis step at the end of the trajectory, the resulting algorithm is exact. Improvements 



of the original update are discussed in [|30|]. 

Here are a few more details: The component of the momentum perpendicular to the X = 
surface is computed using the Feynman-Hellmann theorem SXq/SU = i//q (5/i/5t/)y/b with y/b the 
crossing eigenmode of the kernel operator h and Xq its eigenvalue. Thus 



(3.6) 



TA 



is perpendicular to the surface separating the topological sectors. TA denotes the traceless, anti- 
hermitian component of the vector. With N = A/\A\ the normalized vector perpendicular to the 
surface separating the two topological sectors, the momentum perpendicular to the surface is simply 
7i± = N(n-N). The height of the step AS is computed using Eq. ^j8| which requires one inversion 
of the overlap operator on the crossing mode, independently of how many pseudo-fermion fields 
are used. 
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Figure 1: The history of the topological charge as a function of simulation time. The data is from simulations 
on a 8 4 lattice at lattice spacing a w 0.15 fm [O]. 



How does this algorithm perform at changing topological sector? A time history of the topo- 
logical charge from a simulation on a 8 4 lattice at three values of the bare quark mass is shown in 
Fig. |l|. Whereas the charge changes frequently for the largest quark mass, the simulation tends to 
get stuck in one sector once the quark mass is lowered. As we will see below (Fig. ||) the reason 
is not a lack of attempts to change topology. Rather, the discontinuity is too high to get across. 
It turns out that the distribution of 7Tj_ is not changed from the beginning to the point where the 
trajectory reaches the discontinuity. As shown in Fig. § in the lower right panel, it still has the 
exponential shape exp(— |7r| 2 /2) given by the initial heat bath. This makes values of |7Tj_| 2 above 8 
very unlikely independent of the quark mass. This data, and the one used in the following discus- 
sion comes from a 6 4 lattice with a w 0.15 fm and m 100 MeV. This quark mass is rather large 
and the problems described aggravate once the quark mass is lowered. 

The distribution of the height of the step depends on the quark mass. We expect this because 
a smaller quark mass suppresses configurations of higher topology. However, as it turns out, it 
depends even more on how well the fermion determinant is introduced into the functional integral. 
Fig. ||(top) shows the distribution of the height of the discontinuity changing from v = to |v| = 
1 for one, two and three pseudo-fermion fields. One field yields a very wide distribution and 
most entries have a value above 10 such that tunneling is virtually impossible. Additional pseudo- 
fermions improve the situation; the distribution narrows. That the poor tunneling rate is not a 
problem of physics is shown in the lower left panel of Fig. || where the discontinuity in tr log D^D 
is plotted. If we could compute the fermion determinant and its derivative exactly — without taking 
resort to pseudo-fermions — this is the step we had to overcome. It does not exceed 10 such that 
tunneling is quite likely. The large auto-correlation time in the topological charge is largely (if 
not entirely) due to the poor approximation of the determinant by pseudo-fermions — at least as 
long as the gauge action does not suppress those changes. Fig. || demonstrates that the improved 
estimate indeed helps with the tunneling rate. It shows the fraction of trajectories where the final 
topological charge is different from the initial one as a function of the number of pseudo-fermion 
fields. A roughly linear increase is observed. Note that the cost of a few additional fermions is 
more than recouped by the larger overall step size. That is the reason this improvement method has 
received wide attention in conventional simulations with Wilson type fermions. 

Smaller quark masses make the pseudo-fermion estimate of the fermion determinant even 



7 



Algorithms for dynamical overlap fermions 



Stefan Schaefer 



1 PF 




2 PF 






3 PF 


rlllTTKTh-nni i 




L . . 


r 


i- 





250 500 750 



250 500 750 250 500 750 
2AS 



^1 



-5 5 10 

-2Atr log D + D 




Figure 2: The distribution of the height of the discontinuity in the effective action changing from the v = 
sector to v = ±1 (upper three plots). The lower left plot shows the change in trlogD^D, without recourse 
to pseudo-fermions. The right panel contains the momentum component perpendicular to the surface which 
separates the two topological sectors. The topology is changed if 2A5 > |7r±| 2 . The data is from a HMC 
simulation on 6 4 lattices, a « 0.15fm and am = 0.1. 



more noisy, essentially because the conditioning number of the fermion matrix increases. As the 
previous discussion has shown, this leads to even higher steps and lower tunneling rates. This 
explains why the simulation shown in Fig. |l| — particularly for light quarks — gets stuck for many 
trajectories in sectors of non-zero topology. 

The dependence of the rate of change on the approximation of the determinant might seem 
to interfere with detailed balance. However, detailed balance only makes a statement about the 
ratio of the probability to change from (a configuration in) one topological sector to the other as 
compared to the change in the opposite direction. A better estimate of the determinant, however, 
improves the tunneling rate in both directions and we still get an exact algorithm. The optimal 
tunneling rate is given by the exact determinant which can therefore serve as a guideline on the 
quality of the approximation. We also note that these findings (probably) do not depend on the 
particular algorithm used. Any exact algorithm has the same problems as long as the determinant 
is introduced just by pseudo-fermions. However, as we have shown, different ways to introduce 
the determinant can lead to very different algorithmic performance. 

A particular issue with this algorithm is the scaling with the volume. Each time an eigenmode 
of the kernel operator changes sign, the height of the discontinuity has to be computed. For that the 
overlap operator has to be inverted once on the crossing mode. The cost of this inversion scales at 
least with the volume. Furthermore, the number of attempted crossings turns out to be proportional 
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Figure 3: Left: The dependence of the tunneling rate on the number of pseudo-fermion fields, (data set 
as in Fig. ||) Right: Number of attempted tunnelings per trajectory as a function of the volume. The lattice 
spacing in about a ps 0.15fm, the quark mass m q w lOOMeV. 



to the volume, see Fig. || on the right. This is no surprise since the number of eigenvalues is 
proportional to the volume too. The total cost of the part of the algorithm which deals with changing 
topological sector therefore scales at least with V 2 . 

Whether this scaling law constitutes a problem depends on the actual number of kernel eigen- 
value crossings per trajectory at the target volume. This is mainly a function of the density of 
eigenvalues near the origin which can be reduced by various methods. Gauge actions like the 
Iwasaki action or DBW2 are known to suppress the occurrence of small eigenmodes. This is the 
method of choice in the context of domain wall fermion simulations. It turns out, however, that in 
particular DBW2 has a negative impact on the auto-correlation time of the topological charge [jP]]. 

Another strategy is to use fat links in the definition of the kernel operator. Particularly popular 



are stout links [ 35 ] because they are differentiable and can therefore be used in molecular dynamics 



simulations. So far no impact of the smearing on the number of actual changes of topology has 
been observed. It is believed that the effect of the smearing rather reduces the noise in the motion 
of the eigenmodes and makes it more directed. 

To summarize: The modification of the leap-frog integrator to account for the discontinuity 
due to changing topological sector leads to an exact algorithm. This has been used successfully in 
small volume simulations of QCD. The auto-correlation time of the topological charge is large but 
can be reduced by reducing the noise in the estimator of the fermion determinant. The problematic 
point is that each attempt to change topological sector requires one inversion of the overlap operator. 
The frequency of these attempts scales with the volume which renders this component of the total 
update expensive when the lattice size is increased. 

3.3 Approximate evolution 

In the previous section, we discussed a modification of the leap-frog in the HMC algorithm to 
deal with the step in the action. The discontinuity comes from the sign function in the definition of 
the overlap operator. If instead one uses a continuous approximation to the sign function, the re- 
sulting approximate overlap operator D app is a continuous function of the gauge fields and standard 
HMC can be used. As long as the approximation is close enough to the exact overlap operator, it is 
possible to correct for the difference between D app and D ov at the end. 
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There are two different but related approaches to implement this idea [36]: The first is to 
rewrite the fermion determinant as 

detD 0V = detD a p P detD^D v • (3-7) 

One then simulates the theory with the determinant of D app and reweights at the end to the ex- 
act overlap operator. A variant of this approach is to use the determinant ratio in an additional 
Metropolis step at the end of each trajectory. The second, related possibility is to perform the ini- 
tial heat-bath and the Metropolis step at the end of the trajectory with the exact overlap operator 
whereas the approximation is only used in the computation of the force during the MD evolution. 
To understand the issues involved in this kind of simulation, let us start with the second ap- 



proach put forward by Bode et al. Q37| ] (for a more extensive study see [|36j, [380). They use a fixed 
rational approximation in the definition of the sign function during the evolution. This leads to the 
situation depicted in Fig. |] where the Zolotarov approximation of the sign function is shown such 
that the difference for the part of the spectrum of the kernel operator with |A| > A cut is very small. 
As long as there are not too many eigenmodes with |A| < A cut , the discontinuity in the overlap 
action is approximated by a steep section in 5 e ff which interpolates between the two levels. The 
width of this section is roughly proportional to A cut , the height proportional to AS, the discontinuity 
in the overlap effective action. 

This poses a tuning problem. On the one hand, D app has to be close to the overlap operator 
(Acut has to be small) to easily correct for the difference. On the other hand, a small Ac Ut leads to a 
large derivative of the effective action and to large forces in the molecular dynamics evolution while 
changing topological sector. Therefore a small step-size is required which makes the simulation 



expensive. The analysis in the context of the FKS algorithm of Sec. helps here. The height of 
the discontinuity and thereby the forces in the approximate action can be reduced by improving the 
estimator of the fermion determinant. This leads to smaller forces and an increased probability to 
accept a trajectory during which the topology changed. 

So far the method relies on a stochastic estimate between the determinant of the exact and 
approximate overlap operator. With the specific choice of D app , however, the determinant ratio 



in Eq. \5J] can be computed exactly. It is set up such that the difference between the exact sign 
function and the approximation is zero (or negligible) for |jc| > A cut . From the definition of the 
overlap operator Eq. [2.1[ the following relation is obtained 

N 

£> app = Aw + (*o - T ) £ 8e (A,-) 75 Wi vj ■ (3-8) 

with the sum running over all N eigenmodes \\f of the kernel operator h with an eigenvalue of 
modulus smaller than Ac Ut . The difference between the approximate sign function and the exact 
one is denoted by 8e(x). The ratio between the respective determinant is then given by 

detD a pp = detD 0V det R with R tj = 1 + (R - ^)5e(A,)^— — Wj • (3-9) 

NxN Z Y^i-'ov 

The construction of the NxN matrix R needs N inversions of the overlap operator. This method is 
feasible as long as ,/V is small and the fluctuations in R are moderate. Tests on 8 4 and 10 4 lattices 
are encouraging and subject of an upcoming paper 
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Figure 4: Approximation to the sign function and the effect on the effective action. 



There have been extensive tests using approximate guidance Hamiltonians to simulate the 
overlap operator in the Schwinger model. Christian et al. [[?(]] suggested to use a hypercubic Dirac 
operator in the guidance Hamiltonian. Because the hypercubic operator is already an approximate 
solution of the GW equation, the overlap operator constructed from it is similar to the kernel. 
It turns out that it is not similar enough and the authors find an acceptance rate which vanishes 
with the inverse of the volume. It becomes unacceptably small on their large 32 2 lattices. They 
conclude that a better approximation to the sign function has to be used and discuss possible forms 
in Ref. [38]. Ref. [ |40| ] extends this approach by using in the computation of the force a low order 
polynomial of the hypercubic operator which better approximates the sign function. The volume 
dependence of this approach was not studied. 

Overall, simulating the overlap action using an approximation seems an interesting alterna- 
tive to the modified evolution discussed in the previous section. For differentiable approximations, 
standard methods apply which makes it easier to implement. Also the costly determination of the 
height of the discontinuity where the topology changes is not necessary. However, large forces are 
encountered while changing topological sectors. This requires some tuning to balance ease of sim- 
ulation (small forces) and a small difference between the two operators such that the reweighting is 
possible. How these methods scale with the volume is difficult to answer and needs further study. 
The increased density of modes due to the larger volume probably requires a better approximation 
of the overlap action which in turn is more expensive to simulate. 



3.4 Topology conserving actions 

As discussed in the previous sections, dealing with the discontinuity in the fermionic action 
at the interface of topological sectors is expensive and can lead to a V 2 scaling. An alternative 
approach is to fix the topological charge during the simulation. This eliminates the cost of de- 
termining the height of the step in the FKS algorithm. The fixed topology can be an advantage, 
e.g. in the epsilon regime where predictions are made for fixed topological charge. In particular 
configurations with large charge are virtually impossible to obtain in an algorithm with tunneling 
because they are suppressed by the fermion determinant. Unfortunately, it is unknown whether the 
individual topological sectors are connected. In the following, we will assume that this is the case 
and one can get from any configuration of topology Q to any other by continuous deformation of 
the gauge fields without leaving this particular sector. 
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Still, most physical observables are defined as averages over the full 9 -vacuum. With no ad- 
ditional knowledge about the relative weights of the topological sectors the different sets cannot 
sensibly be combined. Therefore some physics, which is particularly sensitive to the topological 
sampling, cannot be addressed with this method. However, global topology is a boundary condi- 
tion. Thus for most observables, its effect will become negligible for very large volume vanishing 
with 1/V plQ. This is not as rapid as finite volume effects which come from pions winding around 
the torus and exponentially vanish with exp(— m n L). By comparing the results from simulations 
fixed in different topological sectors one can get a handle on these finite volume effects. A re- 
cent quenched study p2|], e.g., found a rather strong dependence of the pseudo-scalar mass on the 
topological charge. 

The obvious way to fix the topology is to eliminate the possibility to tunnel in the FKS al- 
gorithm. One is still left with the reflections but the expensive determination of the height of the 
discontinuity is no longer necessary. An interesting extension of this method [ 13 ] is to compute the 
exact determinant ratios between the two sides of the interface using Eq. ^1]. From this one can 
infer the relative weight of the two topological sectors without need of a stochastic estimate. This 
again requires inversions of the overlap, but the hope is that the relative weights can be measured 
more precisely than from the stochastic estimates by pseudo-fermion fields. 

Another approach which has received considerable attention is to choose an action such that 
the topological sector cannot be changed during molecular dynamics evolution, at least with an 
exact integration of the equations of motion. (The different topological sectors are separated by 
regions of zero weight.) In Ref. [|l4|] it was established that this is the case if all plaquette variables 
are forced to satisfy 

S P [U]<e = ^ (3.10) 



with Sp the plaquette gauge action. Later a looser bound e « could be proved [44]. Unfortu- 
nately, it turned out that the simulation of these gauge actions is difficult [ f45| , |46| ] in particular at 
relevant lattice spacings. 

In a similar spirit is a method suggested by Vranas in the context of domain wall fermions p7|]. 
The idea is to introduce the determinant of the kernel operator det(/i(— Ro)) into the functional inte- 
gral. This determinant is zero where one of the eigenmodes of h(—Ro) changes sign, i.e. where the 
topology as seen by the overlap operator changes. The determinant therefore introduces a repulsive 
potential against the change of the topological sector and eliminates these changes completely if 
one integrates the equations of motion precisely enough. The effect of det(/j(— Ro)) is expected to 
vanish in the continuum limit since it corresponds to a fermion at large negative mass at the cut-off 
scale, essentially acting like a modification of the gauge action. 

JLQCD has started to use a variant of this trick in their large scale simulations with dynamical 
overlap fermions. They use the following ratio of determinants to prevent the trajectories from 
changing topology 

deth 2 (-R ) 



det(/j 2 (-/? () ) +M 2 ) 



(3.11) 



with tunable (twisted) mass parameter }X. This term still is zero only where h(—Ro) has a zero 
mode, however, the effect of the higher modes largely cancels. This cancellation is complete if 
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/i = but then the effect which suppresses the changing of topology vanishes too. The advan- 
tage of fixing the topology in this way is that it proved feasible to simulate (without the overlap 
determinant) in the targeted range of lattice spacings and volumes p8|]. 

At this conference, first results with this method were presented. The simulations were done on 
16 3 x 32 lattices at a lattice spacing around 0. lOfm with various values for the sea quark mass 

4. Summary 

For chiral fermions, the index theorem already holds at finite lattice spacing. This causes the 
fermionic action to be discontinuous where the topological charge defined by the index changes. 
The preferred way to deal with this situation depends on the importance of the topological charge 
in the observable under investigation. Three distinct approaches have been discussed, all based on 
the HMC algorithm. 

The most radical solution is the one discussed last. The action is modified in such a way that 
a continuous modification of the gauge fields does not lead to a change in the topological charge. 
This is achieved by introducing the determinant of the kernel of the overlap operator. Therefore no 
discontinuity is encountered and standard methods can be used to simulate the theory. For most 
observables, the effect of fixed topology is expected to vanish like 1 /V for V —> °°. How large a 
volume is needed for this to be under control will have to be investigated in explicit simulations. 
Such simulations are currently under way and the JLQCD collaboration has already presented first 
results at this conference. 

The other two approaches to dynamical overlap fermions allow for topology to change. The 
first is to modify the molecular dynamics evolution and thereby integrate the discontinuity exactly. 
Each time an eigenmode of the kernel operator changes sign, the height of the step has to be 
determined which costs one inversion of the overlap operator. This leads to V 2 scaling of that part 
of the algorithm because its cost scales with the volume and so does its frequency. The advantage 
is that no large forces occur because of changing topological sector. Also the topological charge 
history is very well under control. 

An alternative approach is to approximate the discontinuity and then correct for this either by 
reweighting or a Metropolis step. This works if one can find an approximate action which on the 
one hand is close enough to the overlap action for the correction term to be small. On the other 
hand, the forces encountered while changing topological sectors have to be small enough to be 
able to integrate the equations of motion to good accuracy with a reasonable step size. This poses 
a tuning problem and again simulations have to show how well this method works in practice. 
Studies in the Schwinger model are encouraging. 

The history of the topological charge is of particular interest in dynamical simulations of chiral 
fermions. We were able to demonstrate that the tunneling rate depends crucially on the way the 
fermion determinant is introduced into the functional integral. Large noise in its estimator leads to 
long auto-correlation times in the topology. So far, only Hasenbusch's mass preconditioning was 
studied and found to be very effective. Other methods might render even better results. 

Dynamical overlap simulations are still very expensive. However, computer power will grow 
and there certainly will be better algorithms too. The benefit of these efforts is a theoretically clean 
description of QCD with chiral symmetry at finite lattice spacing. 
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